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The final black hole left behind after a binary black hole merger can attain a recoil velocity, or a 
“kick”, reaching values up to 5000 km/s. This phenomenon has important implications for gravitational 
wave astronomy, black hole formation scenarios, testing general relativity, and galaxy evolution. 
We consider the gravitational wave signal from the binary black hole merger GW200129_065458 
(henceforth referred to as GW200129), which has been shown to exhibit strong evidence of orbital 
precession. Using numerical relativity surrogate models, we constrain the kick velocity of GW200129 
to vp ~ 15427133, km/s or vf Z 698 km/s (one-sided limit), at 90% credibility. This marks the 
first identification of a large kick velocity for an individual gravitational wave event. Given the 
kick velocity of GW200129, we estimate that there is a less than 0.48% (7.7%) probability that the 
remnant black hole after the merger would be retained by globular (nuclear star) clusters. Finally, 
we show that kick effects are not expected to cause biases in ringdown tests of general relativity for 
this event, although this may change in the future with improved detectors. 


Introduction.— When two black holes (BHs) orbit each 
other, they emit gravitational waves (GWs) which carry 
away energy and angular momentum. This causes the 
orbit to shrink in a runaway process that culminates in 
the merger of the BHs into a single remnant BH. At the 
same time, the GWs can also carry away linear momentum 
from the binary, shifting its center of mass in the opposite 
direction [1]. Most of the linear momentum is lost near 
the merger [2], resulting in a recoil or “kick” velocity 
imparted to the remnant BH. 

Kicks are particularly striking for precessing binaries, 
in which the component BH spins are tilted with respect 
to the orbital angular momentum. For these systems, the 
spins interact with the orbital angular momentum as well 
as with each other, causing the orbital plane to precess [3]. 
Numerical relativity (NR) simulations revealed that the 
kick velocities for precessing binaries can reach values up 
to ~ 5000 km/s [4-6], large enough to be ejected from 
any host galaxy [7]. 


Kicks have important implications for BH astrophysics. 


Following a supermassive BH merger, the remnant BH can 
be displaced from the galactic center or ejected entirely [7], 
impacting the galaxy’s evolution [8], fraction of galaxies 
with central supermassive BHs [9], and event rates [10] 
for the future LISA mission [11]. For stellar-mass BHs 
like those observed by LIGO [12] and Virgo [13], kicks can 
limit the formation of heavy BHs. BH masses greater than 
~65Mo are disfavored by supernova simulations [14, 15], 
but have been seen in GW events [16-18]. This could be 
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explained by second-generation mergers [19], in which one 
of the component BHs is itself a remnant from a previous 
merger, and is thus more massive than the original stellar- 
mass progenitors. However, if the kick from the first 
merger is large enough, the remnant BH would get ejected 
from its host galaxy and would not participate in another 
merger. 


Unfortunately, observational evidence of large kicks has 
been elusive. While various candidates from electromag- 
netic observations have been identified, their nature is 
debated [20]. Similarly, observing kicks using GW sig- 
nals has been challenging [21-27]. For example, Varma 
et al. [21] used accurate models based on NR simulations 
to show that kicks from precessing binaries can be reli- 
ably inferred with LIGO-Virgo operating at their design 
sensitivity. However, the GW events analyzed in Ref. [21], 
which only included signals in the first two LIGO-Virgo 
observing runs [28], were not loud enough to constrain 
the kick. 


Since then, the LIGO-Virgo detectors have been further 
upgraded, and the GW data from the third observing 
run were released in two stages, O8a [17] and O3b [18]. 
Notably, O3a provided the first evidence for precession 
in the ensemble population of merging binaries [29], even 
though none of the individual GW events unambiguously 
exhibited precession [17]. Finally, in O3b, the binary BH 
merger GW200129 was identified as the first individual 
GW event showing strong evidence of precession [18, 30]. 
Similarly, support for large kicks was identified in the 
ensemble population using the O8a data [31, 32], even 
though the individual events were not loud enough for an 
unambiguous kick inference [21, 33] (with the exception 


of GW190814 [34], which was found to have a small kick 
of ~ 74+7° km/s at 90% credibility [27]). 

In this Letter, we use the method developed in 
Ref. [21] to show that GW200129 has a large kick ve- 
locity (~ 1542*755g km/s at 90% credibility). As an ap- 
plication of the kick constraint, we compute the retention 
probability for the remnant BH of GW200129 in various 
host environments, and discuss the implications for the 
formation of heavy stellar-mass BHs. Finally, we show 
that Doppler effects due to the kick on the remnant mass 
measurement are small for this event, and should not 
impact ringdown tests of general relativity (GR). 


Methods.— We follow the procedure outlined in Ref. [21] 
to infer the kick from a GW signal. We begin by measur- 
ing the binary source parameters following Bayes’ theo- 
rem [35]: 


(Ad) x L(d|A) (A), (1) 


where p(A|d) is the posterior probability distribution of 
the binary parameters A given the observed data d, £(d|X) 
is the likelihood of the data given A, and 7(A) is the prior 
probability distribution for A. Under the assumption of 
Gaussian detector noise, the likelihood £(d|X) can be eval- 
uated for any A using a gravitational waveform model and 
the observed data stream d [35]. A stochastic sampling 
algorithm is then used to draw posterior samples for A 
from p(A|d). We use the Parallel Bilby [36] parameter 
estimation package with the dynesty [37] sampler. 

For quasicircular binary BHs, the full set of parameters 
A is 15 dimensional [18]. This includes the 8D intrinsic 
parameters: the component masses (mı and m2) and 
spins (xı and x2, each of which is a 3D vector), as 
well as the 7D extrinsic parameters: the distance, right 
ascension, declination, time of arrival, coalescence phase, 
binary inclination, and polarization angle. Here, index 
1 (2) corresponds to the heavier (lighter) BH, x1,2 are 
dimensionless spins with magnitudes x 1,2 < 1, and masses 
refer to the detector frame redshifted masses. We also 
define the mass ratio q = m1/m2 > 1, total mass M = 
mı + mg, and use geometric units with G = c= 1. 

We employ the NR surrogate models NRSur7dq4 [38] 
and NRSur7dq4Remnant [38, 39] to infer the kick. Con- 
structed by effectively interpolating between ~ 1500 pre- 
cessing NR simulations, NRSur7dq4 predicts the gravita- 
tional waveform, while NRSur7dq4Remnant predicts the 
mass my, spin x, and kick velocity vy of the remnant BH. 
We first obtain posterior samples for all 15 binary parame- 
ters using NRSur7dq4. The spins are measured in a source 
frame defined at a given reference point (see below): the 
z-axis lies along the instantaneous orbital angular momen- 
tum, the z-axis points along the line of separation from 
the lighter to the heavier BH, and the y-axis completes the 
right-handed triad. The remnant properties, which are 
also defined in the same source frame, depend only on the 
intrinsic parameters A = {m1, M2, X1; X2}. Therefore, 
the posteriors for the remnant properties are obtained by 
evaluating NRSur7dq4Remnant on the A posterior samples; 


put simply, NRSur7dq4Remnant is a function of A that 
yields my, xs and vy. We also compute the effective prior 
distribution for vy, by evaluating NRSur7dq4Remnant on 
A samples drawn from the prior r(A). The difference 
between the kick posterior and prior can be used to gauge 
how informative the data are about the kick [21]. 


Traditional modeling methods assume a phenomenolog- 
ical ansatz for the waveform [40, 41] or remnant proper- 
ties [42-45], and calibrate remaining free parameters to 
NR simulations. NR surrogate methods [38, 39, 46, 47], 
on the other hand, take a data-driven approach and the 
models are trained directly against precessing NR simu- 
lations. In this approach, one first constructs a suitable 
numerical basis using a subset of the NR waveforms, 
and then builds fits across parameter space for the ba- 
sis coefficients; we refer the reader to Ref. [38] for more 
details. NR surrogate models do not need to introduce 
additional assumptions about the underlying phenomenol- 
ogy which would necessarily introduce some systematic 
error. Through cross-validation studies, it has been shown 
that both NRSur7dq4 and NRSur7dq4Remnant achieve ac- 
curacies comparable to the simulations themselves [38], 
and as a result, are the most accurate models currently 
available for precessing systems, within their parameter 
space of validity: both models are trained on simulations 
with q < 4 and %1,2 < 0.8, but can be extrapolated to 
q <6 and %1,2 < 1 [38]. As GW200129 shows significant 
support for large spins, we conduct some tests of the 
surrogate models in this regime in the Supplement [48]. 


For the prior in Eq. (1), we follow Ref. [18] and adopt 
a uniform prior for spin magnitudes (with 0 < y1, x2 < 
0.99) and redshifted component masses, an isotropic prior 
for spin orientations, sky location and binary orientation, 
and a distance prior (UniformSourceFrame [49]) that as- 
sumes uniform source distribution in comoving volume 
and time. In addition, we place the following constraints: 
q < 6 and 60 < M < 400. These constraints are moti- 
vated by the regime of validity of NRSur7dq4 [38], and are 
broad enough to safely encompass the posterior spread of 
GW 200129 [18]. 

Because the spin directions are not constant for pre- 
cessing binaries, spin measurements are inherently tied 
to a specific moment in the binary’s evolution. The stan- 
dard approach is to measure the spins at the point where 
the frequency of the GW signal at the detector reaches 
a prespecified reference value, typically fret =20 Hz [17]. 
This is mainly motivated by the fact that the sensitivity 
band of current detectors begins near this value [12, 13]. 
However, Ref. [50] recently showed that constraints on 
orbital-plane spin directions can be greatly improved by 
measuring the spins near the merger, in particular, at a 
fixed dimensionless reference time tret /M =—100 before 
the peak of the GW amplitude. This improvement can 
be attributed to the waveform being more sensitive to 
variations in the orbital-plane spin directions near the 
merger [50] (tre/M =—100 typically falls within ~ 2 — 4 
GW cycles before the peak amplitude, independent of the 
binary parameters [50].). 
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Figure 1. Constraints on the mass ratio and spins for 
GW 200129, at reference time trer/M =—100. The dark (light) 
regions represent the 50% (90%) credible bounds on joint 2D 
posteriors, while the diagonal plots show 1D marginalized pos- 
teriors. There is a preference for large xı and cos 4; ~ 0, mean- 
ing there is substantial spin in the orbital plane, which leads to 
precession. For comparison, we also show the 1D marginalized 
posteriors at fref = 20 Hz. The azimuthal spin angles (espe- 
cially 1) are much better constrained at tref /M = —100; this 
is critical for constraining the kick. 


We will adopt the tref /M = —100 reference point for the 
main results in this paper, but will show a comparison 
against fret = 20 Hz for completeness. As we will dis- 
cuss below, the choice of reference point has a negligible 
impact on the kick inference itself, but comparing the 
spin posteriors at the two reference points helps illustrate 
why a kick constraint is possible in the first place. As 
a bonus, spin measurements at tye¢/M =—100 are con- 
venient for inferring the kick as the NRSur7dq4Remnant 
model is also trained at this reference time [38]; this choice 
was found to lead to a more accurate remnant BH model 
in Refs. [38, 39]. 


GW200129 spin measurements.— GW200129 is the 
first GW event showing strong signs of precession [18, 30]. 
Figure 1 shows the posterior distribution for the mass 
ratio and spin parameters obtained using the NRSur7dq4 
model at reference points tref /M =—100 and fref = 20 Hz; 
our constraints at fret =20 Hz are consistent with those 
of Ref. [30]. The spin vectors x1,2 are decomposed into 
magnitudes x12, tilts angles 01,2 with respect to the z-axis, 
and azimuthal angles ¢1,2 with respect to the x-axis of 
the source frame. Due to precession, spins measurements 
vary between the two reference points but can be related 
by a spin evolution [38, 50]. 


For both reference points in Fig. 1, there is a clear 
preference for large orbital-plane spins for the heavier 
BH (large xı and cos; ~ 0). Even though the spin of 
the lighter BH is not well measured, this is sufficient for 
precession. We stress that while precessing binaries tend 
to have larger kicks [4-6], precession does not necessarily 
imply a large kick, and it is important to directly compute 
the kick velocity as we do in the next section. In particular, 
the kick can vary from zero to ~ 5000 km/s just by 
changing the azimuthal spin angles, even for systems with 
large orbital-plane spins [4—6]. 

Next, the azimuthal angles (especially ¢,) in Fig. 1 
are much better constrained at tyer/M =—100, while the 
other parameters do not change significantly. | This fea- 
ture is key: even though the azimuthal angles are poorly 
constrained in the inspiral, they are well constrained at 
tre¢/M = —100 [50]. As the kick depends sensitively on 
the azimuthal angles near the merger [4], successfully 
measuring these angles at tre¢/M = —100 is critical for 
constraining the kick. 

Spins measured at fret = 20 Hz can also be evolved 
consistently to tres/M = —100 using NRSur7dq4 dy- 
namics [38]. In fact, this procedure is internally ap- 
plied by NRSur7dq4Remnant if the spins are specified at 
fret =20 Hz [38, 39]. Therefore, by construction, the kick 
posterior for individual GW events is independent of the 
reference point at which the spins are initially measured 
(modulo NRSur7dq4 spin evolution errors, which are small 
compared to the model errors [38, 47]). Therefore, for 
the purpose of this paper, the main benefit of the spin 
measurements at tre¢/M = —100 is to illustrate why a 
successful kick constraint is possible in the first place. 
The supplement of Ref. [31] discusses other benefits, in 
particular, for constraining the ensemble population of 
spins and kicks. In the rest of the paper, we will use the 
spin measurements at tref / M =—100. 

As noted in Refs. [18, 30], the inference of preces- 
sion in GW200129 depends on the waveform model 
used. In particular, while the phenomenological model 
IMRPhenomXPHM [40] recovers precession, the effective-one- 
body model SEOBNRV4PHM [41] does not. Among these 
models, only NRSur7dq4 is informed by precessing NR sim- 
ulations and is more accurate by about an order of magni- 
tude [38]. By contrast, SEOBNRv4PHM and IMRPhenomXPHM 
approximate precession effects by “twisting” the frame 
of an equivalent aligned-spin binary [40, 41]. Further- 
more, Ref. [50] found that NRSur7dq4 is necessary to 
accurately measure the spin vectors x1,2, in particular, 
the spin directions within the orbital plane [50], which 
have a strong influence on the kick [4]. Similarly, given 
the spin measurements, NRSur7dq4Remnant is necessary 
to accurately predict the kick velocity [21]. For these 


1 The posteriors for x1, X2, and q are expected to be consistent 
between the two reference points (modulo parameter estimation 
uncertainty) as these parameters are independent of the reference 
point. 
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Figure 2. Left: Kick magnitude constraints for GW200129. We show the posterior and the effective prior, along with known 
ranges for the escape velocities for various types of host environments for comparison. There is a clear preference for large kicks 
in the posterior, with vf 2 698 km/s at 90% credibility. Right: Cumulative distribution functions (CDFs) for the kick posterior 
and prior. The upper bounds of the escape velocity ranges from the left panel are shown as vertical dotted lines. The upper 
limit for retention probability of the merger remnant is given by the intersection of these lines with the posterior CDF. 


reasons, we treat NRSur7dq4 and NRSur7dq4Remnant as 
the preferred models for analyzing GW200129. 


GW200129 kick velocity.— Figure 2 shows our con- 
straints on the kick magnitude vf of GW200129, obtained 
by evaluating NRSur7dq4Remnant on the NRSur7dq4 A 
posteriors at trep/M =—100. In the left panel, we show 
the posterior and prior distributions for vy, along with 
fiducial escape velocities for globular clusters [51], nu- 
clear star clusters [51], giant elliptical galaxies [7] and 
Milky Way-like galaxies [52] for comparison. Unlike the 
events considered in Ref. [21], the vy posterior is clearly 
distinguishable from the prior, and there is substantial 
information gain about the kick.” 

The kick magnitude is constrained to vf ~ 1543ta 
km/s (median and 90% symmetric credible interval), or 
vf Z 698 km/s (lower 10th percentile), making GW200129 
the first GW event identified as having a large kick velocity. 
We note, however, that such large kick velocities are not 
surprising given previous constraints on the ensemble 
properties of merging binary BHs [31, 32]. For example, 
Fig. 3 of Ref. [31], which shows estimates of the ensemble 
kick distribution, includes nonneglibile support up to 
vf ~ 1500 km/s. 

The large kick of GW200129 raises the question of 
whether the remnant BH is ejected from its host environ- 
ment. This has implications for the formation of heavy 


2 The Kullback-Leibler (KL) divergence [53] from the prior to 
the posterior in Fig. 2 is 4.3 bits. By contrast, the largest KL 
divergence for the events considered in Ref. [21] was 0.22 bits. 


BHs through second-generation mergers in dense environ- 
ments [19]. This formation channel is one possible way 
to explain observations of BHs with masses = 65Mọ [16- 
18], which fall within the mass gap expected due to the 
(pulsational) pair-instability supernova processes [14, 15]. 
To address this, we compute the retention probability for 
the remnant BH of GW200129 in globular clusters and 
nuclear star clusters, both of which host dense stellar envi- 
ronments where merger remnants can potentially interact 
with other BHs and form binaries. 


The right panel of Fig. 2 shows the cumulative distribu- 
tion functions (CDFs) for the vy posterior and prior. As 
the posterior CDF (vy) denotes the probability that the 
kick magnitude of GW200129 is below vy, we take it to be 
the probability that the remnant BH is retained by a host 
environment with an escape velocity of vy. The vertical 


dotted lines indicate the maximum escape velocity v232* 


for various host environments; CDF(v%2*) sets the upper 
limit on the retention probability for that host. In partic- 
ular, assuming v%*=100 km/s (v22* =600) for globular 
(nuclear star) clusters there is a less than 0.48% (7.7%) 
probability that the remnant BH of GW200129 is retained 
by those hosts. This is consistent with Refs. [31,32], where 
globular clusters were already identified as an unlikely site 
for second generation mergers, even for more moderate 


kicks. 


Remnant mass and Doppler shifts.— Our method 
provides predictions for both the magnitude and direction 
of the kick [21]. If the kick vector vy has a significant com- 
ponent along (or opposite) the line-of-sight, the observed 
GW signal can be influenced by the kick. At leading order, 


00 02 04 06 08 1.0 00 02 04 06 08 1.0 
Normalized probability density Normalized probability density 


Figure 3. Posterior samples for the full kick vector vy in the 
source frame at tret /M =—100. Each purple marker indicates 
a kick posterior sample; an arrow drawn from the origin to the 
marker would show the kick vector vs. The outer radius of the 
sphere corresponds to vf = 2500 km/s. The z-axis (orange) 
and y-axis (green) are shown as arrows near the origin; the x—y 
plane is orthogonal to the orbital angular momentum direction. 
The blue markers on the sphere show posterior samples for the 
line-of-sight direction to the observer. For both distributions, 
the spread represents the measurement uncertainty, and the 
color reflects posterior probability density (normalized so that 
the peak density is 1). A rotating perspective of this plot can 
be seen at vijayvarma392.github.io/GW200129/#kick. 


the kick’s effect can be described as a Doppler shift of the 
GW frequency [23]. However, as GR lacks any intrinsic 
length scales, a uniform increase in signal frequency is 
completely degenerate with a decrease in total mass M, 
and vice versa. Thus, if not explicitly accounted for, a 
frequency shift due to a kick can bias mass measurements. 
In particular, because the kick is mostly imparted near 
the merger [2], the Doppler shift only affects the merger 
and ringdown part of the signal. This can lead to biases 
in the measurement of the remnant mass my, [21, 54], 
and potentially impact tests of GR using the ringdown 
signal [55]. However, this effect is expected to be small 
for current detectors [21, 23]. 

In the following, we verify that the Doppler effect on 
the remnant mass measurement of GW200129 is indeed 
small. At leading order, the Doppler-shifted remnant 
mass is given by [23]: 


mP’ = my (l+vy-/c), (2) 


where c is the speed of light and f is the unit vector 
pointing along the line-of-sight from the observer to the 
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Figure 4. The remnant mass and the Doppler shifted remnant 
mass for GW200129, as inferred in the detector frame. There 
is an overall redshift, as the kick direction in Fig. 3 is pointed 
(roughly) away from the observer. However, as these distribu- 
tions are very close, we do not expect ringdown tests of GR 
to be impacted by the kick for this event. 


source. The line-of-sight direction is obtained from our 
inference setup, parameterized by (+, ¢). ¢ is the inclina- 
tion angle between the orbital angular momentum and 
the line-of-sight to the observer, and ¢ is the azimuthal 
angle to the observer in the orbital plane, both defined in 
the source frame at tret /M =—100. 

Figure 3 shows the posterior distributions for the 
full kick vector vy (also defined in the source frame at 
tret/M = —100) and the line-of-sight direction. We find 
that the kick and the line-of-sight are not very well (anti-) 
aligned; therefore, we do not expect significant Doppler 
shifts for this signal. Finally, Fig. 4 shows the posterior 
distributions for m, (obtained from NRSur7dq4Remnant) 
and mS (computed using Eq. (2)) for GW200129. As ex- 
pected, the difference between these distributions is very 
small compared to the measurement uncertainty, meaning 
that tests of GR should not be impacted by the Doppler 
effect for this event. As detector sensitivity improves, this 
may not be the case, however, and it may be necessary 
to explicitly account for this effect [21]. 


Conclusions.— We use NR surrogate models for the 
gravitational waveform and the remnant BH proper- 
ties to infer the kick velocity for the binary BH merger 
GW200129. The kick magnitude is constrained to vf ~ 
1542+767, km/s or vf > 698 km/s, at 90% credibility. 
Given the kick velocity, we estimate that there is at 
most a 0.48% (7.7%) probability that the remnant BH of 
GW200129 would be retained by globular (nuclear star) 
clusters. Finally, we show that the Doppler effect on the 
remnant mass is small compared to current measurement 
uncertainty; therefore ringdown tests of GR are not ex- 
pected to be significantly impacted by the kick for this 
event. 

Observational evidence for kicks has far reaching im- 
plications for BH astrophysics. GW200129 is the first 


GW event identified as having a large kick velocity. Large 
kicks like this have been previously predicted based on the 
ensemble kick distribution of merging binary BHs [31, 32], 
and we can expect to see more such events as detector 
sensitivity improves. In particular, such observations 
can help resolve the mystery of the heavy BHs seen by 
LIGO-Virgo [16-18], by constraining the rate of second- 
generation mergers. 


Acknowledgments.— We thank Arif Shaikh for com- 
ments on the manuscript. V.V acknowledges funding 
from the European Union’s Horizon 2020 research and 
innovation program under the Marie Skłodowska-Curie 
grant agreement No. 896869. S.B., C.-J.H., and S.V. 
acknowledge support of the National Science Founda- 
tion (NSF) and the LIGO Laboratory. S.B. is also sup- 
ported by the NSF Graduate Research Fellowship un- 
der Grant No. DGE-1122374. S.V. is also supported by 


[1] M. J. Fitchett, “The 
wave momentum losses 


influence of gravitational 

on the centre of mass 

motion of a Newtonian binary system,” Monthly 

Notices of the Royal Astronomical Society 203, 

1049-1062 (1983), http://oup.prod.sis.lan/mnras/article- 

pdf/203/4/1049/18223796/mnras203-1049.pdf. 

Jose A. Gonzalez, Ulrich Sperhake, Bernd Bruegmann, 

Mark Hannam, and Sascha Husa, “Total recoil: The Max- 

imum kick from nonspinning black-hole binary inspiral,” 

Phys. Rev. Lett. 98, 091101 (2007), arXiv:gr-qc/0610154 

[gr-qc]. 

Theocharis A. Apostolatos, Curt Cutler, Gerald J. Suss- 

man, and Kip S. Thorne, “Spin-induced orbital precession 

and its modulation of the gravitational waveforms from 

merging binaries,” Phys. Rev. D 49, 6274-6297 (1994). 

[4] Manuela Campanelli, Carlos O. Lousto, Yosef Zlochower, 
and David Merritt, “Maximum gravitational recoil,” Phys. 
Rev. Lett. 98, 231102 (2007), arXiv:gr-qc/0702133 [GR- 
QC. 

[5] J. A. Gonzalez, M. D. Hannam, U. Sperhake, Bernd 
Bruegmann, and S. Husa, “Supermassive recoil veloci- 
ties for binary black-hole mergers with antialigned spins,” 
Phys. Rev. Lett. 98, 231101 (2007), arXiv:gr-qc/0702052 
[GR-QC]. 

[6] Carlos O. Lousto and Yosef Zlochower, “Hangup Kicks: 
Still Larger Recoils by Partial Spin/Orbit Alignment of 
Black-Hole Binaries,” Phys. Rev. Lett. 107, 231102 (2011), 
arXiv:1108.2009 [gr-qc]. 

[7] David Merritt, Milos Milosavljevic, Marc Favata, Scott A. 
Hughes, and Daniel E. Holz, “Consequences of gravita- 
tional radiation recoil,” Astrophys. J. 607, L9-L12 (2004), 
ar Xiv:astro-ph/0402057 [astro-ph]. 

[8] S. Komossa and David Merritt, “Gravitational Wave Re- 

coil Oscillations of Black Holes: Implications for Unified 

Models of Active Galactic Nuclei,” Astrophys. J. 689, 

L89 (2008), arXiv:0811.1037 [astro-ph]. 

Marta Volonteri, Kayhan Gültekin, and Massimo Dotti, 

“Gravitational recoil: effects on massive black hole occu- 

pation fraction over cosmic time,” Monthly Notices of 


[2 


[3 


(9 


NSF Grant No. PHY-2045740. T.I. is supported by the 
Heising-Simons Foundation, the Simons Foundation, and 
NSF Grants Nos. PHY-1748958, PHY-1806665 and DMS- 
1912716. F.S. and S.E.F. are supported by NSF Grants 
Nos. PHY-2110496 and PHY-1806665. Computations 
were performed on the Wheeler cluster at Caltech, which 
is supported by the Sherman Fairchild Foundation and 
by Caltech; and the High Performance Cluster at Caltech. 
This material is based upon work supported by NSF’s 
LIGO Laboratory which is a major facility fully funded 
by the NSF. LIGO was constructed by the California 
Institute of Technology and Massachusetts Institute of 
Technology with funding from the NSF and operates un- 
der cooperative agreement PHY-0757058. This research 
made use of data, software and/or web tools obtained 
from the Gravitational Wave Open Science Center [56], 
a service of the LIGO Laboratory, the LIGO Scientific 
Collaboration and the Virgo Collaboration. 


the Royal Astronomical Society 404, 2143-2150 (2010), 
arXiv:1001.1743 [astro-ph.CO]. 

10] A. Sesana, “Extreme recoils: impact on the detec- 

tion of gravitational waves from massive black hole bi- 

naries,” Mon. Not. Roy. Astron. Soc. 382, 6 (2007), 

arXiv:0707.4677 [astro-ph]. 

11] Pau Amaro-Seoane et al., “Laser Interferometer Space 

Antenna,” arXiv e-prints , arXiv:1702.00786 (2017), 

ar Xiv:1702.00786 [astro-ph.IM]. 

12] J. Aasi et al. (LIGO Scientific), “Advanced LIGO,” Class. 

Quant. Grav. 32, 074001 (2015), arXiv:1411.4547 [gr-qc]. 

13] F. Acernese et al. (Virgo), “Advanced Virgo: a second- 

generation interferometric gravitational wave detector,” 

Class. Quant. Grav. 32, 024001 (2015), arXiv:1408.3978 

[gr-qc]. 

14] S. E. Woosley, “Pulsational Pair-Instability Supernovae,” 

Astrophys. J. 836, 244 (2017), arXiv:1608.08939 [astro- 

ph.HE}. 

15] Pablo Marchant, Mathieu Renzo, Robert Farmer, Kaliroe 
M. W. Pappas, Ronald E. Taam, Selma de Mink, and Vas- 
siliki Kalogera, “Pulsational pair-instability supernovae in 

very close binaries,” (2018), 10.3847/1538-4357 /ab3426, 
arXiv:1810.13412 [astro-ph.HE]. 

[16] R. Abbott et al. (LIGO Scientific, Virgo), “GW190521: A 
Binary Black Hole Merger with a Total Mass of 150 Mo,” 

Phys. Rev. Lett. 125, 101102 (2020), arXiv:2009.01075 
[gr-qc]. 

[17] R. Abbott et al. (LIGO Scientific, Virgo), “GWTC-2: 
Compact Binary Coalescences Observed by LIGO and 
Virgo During the First Half of the Third Observing Run,” 
Phys. Rev. X 11, 021053 (2021), arXiv:2010.14527 [gr-qc]. 

[18] R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), 
“GWTC-3: Compact Binary Coalescences Observed by 
LIGO and Virgo During the Second Part of the Third 

Observing Run,” (2021), arXiv:2111.03606 [gr-qc]. 

[19] Davide Gerosa and Maya Fishbach, “Hierarchical mergers 
of stellar-mass black holes and their gravitational-wave 
signatures,” Nature Astron. 5, 8 (2021), arXiv:2105.03439 
[astro-ph. HE]. 


[20] 


[21] 


[22] 


[23] 


[24] 


[25] 


[26] 


[27] 


[28] 


[29] 


[30] 


[31] 


[32] 


[33] 


[34] 


Tamara Bogdanovic, M. Coleman Miller, and Laura 
Blecha, “Electromagnetic Counterparts to Massive Black 
Hole Mergers,” (2021), arXiv:2109.03262 [astro-ph.HE]. 
Vijay Varma, Maximiliano Isi, and Sylvia Biscoveanu, 
“Extracting the Gravitational Recoil from Black Hole 
Merger Signals,” Phys. Rev. Lett. 124, 101104 (2020), 
ar Xiv:2002.00296 [gr-qc]. 

R. Abbott et al. (LIGO Scientific, Virgo), “Properties and 
Astrophysical Implications of the 150 Mo Binary Black 
Hole Merger GW190521,” Astrophys. J. 900, L13 (2020), 
ar Xiv:2009.01190 [astro-ph.HE]. 

Davide Gerosa and Christopher J. Moore, “Black hole 
kicks as new gravitational wave observables,” Phys. Rev. 
Lett. 117, 011101 (2016), arXiv:1606.04226 [gr-qc]. 
Juan Calderén Bustillo, James A. Clark, Pablo Laguna, 
and Deirdre Shoemaker, “Tracking black hole kicks from 
gravitational wave observations,” Phys. Rev. Lett. 121, 
191102 (2018), arXiv:1806.11160 [gr-qc]. 

Carlos O. Lousto and James Healy, “Kicking gravitational 
wave detectors with recoiling black holes,” Phys. Rev. 
D100, 104039 (2019), arXiv:1908.04382 [gr-qc]. 

James Healy, Carlos O. Lousto, Jacob Lange, and Richard 
O’Shaughnessy, “Application of the third RIT binary 
black hole simulations catalog to parameter estimation of 
gravitational waves signals from the LIGO-Virgo 01/02 
observational runs,” Phys. Rev. D 102, 124053 (2020), 
ar Xiv:2010.00108 [gr-qc]. 

Parthapratim Mahapatra, Anuradha Gupta, Marc Favata, 
K. G. Arun, and B. S. Sathyaprakash, “Remnant Black 
Hole Kicks and Implications for Hierarchical Mergers,” 
Astrophys. J. Lett. 918, L31 (2021), arXiv:2106.07179 
[astro-ph. HE]. 

B. P. Abbott et al. (LIGO Scientific, Virgo), “GWTC- 
1: A Gravitational-Wave Transient Catalog of Compact 
Binary Mergers Observed by LIGO and Virgo during 
the First and Second Observing Runs,” Phys. Rev. X9, 
031040 (2019), arXiv:1811.12907 [astro-ph.HE]. 

R. Abbott et al. (LIGO Scientific, Virgo), “Population 
Properties of Compact Objects from the Second LIGO- 
Virgo Gravitational-Wave Transient Catalog,” Astrophys. 
J. Lett. 913, L7 (2021), arXiv:2010.14533 [astro-ph.HE]. 
Mark Hannam, Charlie Hoy, Jonathan E. Thompson, 
Stephen Fairhurst, Vivien Raymond, and members of the 
LIGO (VIRGO), “Measurement of general-relativistic pre- 
cession in a black-hole binary,” (2021), arXiv:2112.11300 
[gr-qc]. 

Vijay Varma, Sylvia Biscoveanu, Maximiliano Isi, Will M. 
Farr, and Salvatore Vitale, “Hints of spin-orbit reso- 
nances in the binary black hole population,” (2021), 
ar Xiv:2107.09693 [astro-ph.HE]. 

Zoheyr Doctor, Ben Farr, and Daniel E. Holz, “Black 
Hole Leftovers: The Remnant Population from Binary 
Black Hole Mergers,” Astrophys. J. Lett. 914, L18 (2021), 
ar Xiv:2103.04001 [astro-ph.HE]. 

Tousif Islam, Feroz Shaik, Carl-Johan Haster, Vijay 
Varma, Scott Field, Jacob Lange, Richard O’Shaughnessy, 
Rory Smith, and Avi Vajpeyi, “Re-analysis of GWTC- 
3 events with Numerical Relativity Surrogate models,” 
(2022), in preparation. 

R. Abbott et al. (LIGO Scientific, Virgo), “GW190814: 
Gravitational Waves from the Coalescence of a 23 Solar 
Mass Black Hole with a 2.6 Solar Mass Compact Object,” 
Astrophys. J. Lett. 896, L44 (2020), arXiv:2006.12611 
[astro-ph. HE]. 


[35] 


[36] 


[37] 


[38] 


[39] 


[40] 


[41] 


[42] 


[43] 


[44] 


[45] 


[46] 


[47] 


[48] 


Eric Thrane and Colm Talbot, “An introduction to 
Bayesian inference in gravitational-wave astronomy: Pa- 
rameter estimation, model selection, and hierarchical mod- 
els,” Publications of the Astronomical Society of Australia 
36, e010 (2019), arXiv:1809.02293 [astro-ph.IM]. 

Rory J.E. Smith, Gregory Ashton, Avi Vajpeyi, and 
Colm Talbot, “Massively parallel Bayesian inference for 
transient gravitational-wave astronomy,” Mon. Not. Roy. 
Astron. Soc. 498, 4492-4502 (2020), arXiv:1909.11873 
[gr-qc]. 

Joshua S. Speagle, “DYNESTY: a dynamic nested sam- 
pling package for estimating Bayesian posteriors and ev- 
idences,” Monthly Notices of the Royal Astronomical 
Society 493, 3132-3158 (2020), arXiv:1904.02180 [astro- 
ph.IM]. 

Vijay Varma, Scott E. Field, Mark A. Scheel, Jonathan 
Blackman, Davide Gerosa, Leo C. Stein, Lawrence E. Kid- 
der, and Harald P. Pfeiffer, “Surrogate models for precess- 
ing binary black hole simulations with unequal masses,” 
Phys. Rev. Research. 1, 033015 (2019), arXiv:1905.09300 
[gr-qc]. 

Vijay Varma, Davide Gerosa, Leo C. Stein, François 
Hébert, and Hao Zhang, “High-accuracy mass, spin, and 
recoil predictions of generic black-hole merger remnants,” 
Phys. Rev. Lett. 122, 011101 (2019), arXiv:1809.09125 
[gr-qc]. 

Geraint Pratten et al., “Computationally efficient models 
for the dominant and subdominant harmonic modes of 
precessing binary black holes,” Phys. Rev. D 103, 104056 
(2021), arXiv:2004.06503 [gr-qc]. 

Serguei Ossokine et al., “Multipolar Effective-One-Body 
Waveforms for Precessing Binary Black Holes: Construc- 
tion and Validation,” Phys. Rev. D 102, 044055 (2020), 
arXiv:2004.09442 [gr-qc]. 

Fabian Hofmann, Enrico Barausse, and Luciano Rezzolla, 
“The final spin from binary black holes in quasi-circular 
orbits,” Astrophys. J. 825, L19 (2016), arXiv:1605.01938 
[gr-qc]. 

Enrico Barausse, Viktoriya Morozova, and Luciano Rez- 
zolla, “On the mass radiated by coalescing black-hole 
binaries,” Astrophys. J. 758, 63 (2012), [Erratum: Astro- 
phys. J.786,76(2014)], arXiv:1206.3803 [gr-qc]. 

Xisco Jiménez-Forteza, David Keitel, Sascha Husa, Mark 
Hannam, Sebastian Khan, and Michael Piirrer, “Hierar- 
chical data-driven approach to fitting numerical relativity 
data for nonprecessing binary black holes with an applica- 
tion to final spin and radiated energy,” Phys. Rev. D95, 
064024 (2017), arXiv:1611.00332 [gr-qc]. 

Carlos O. Lousto, Yosef Zlochower, Massimo Dotti, and 
Marta Volonteri, “Gravitational Recoil From Accretion- 
Aligned Black-Hole Binaries,” Phys. Rev. D85, 084015 
(2012), arXiv:1201.1923 [gr-qc]. 

Vijay Varma, Scott E. Field, Mark A. Scheel, Jonathan 
Blackman, Lawrence E. Kidder, and Harald P. Pfeiffer, 
“Surrogate model of hybridized numerical relativity binary 
black hole waveforms,” Phys. Rev. D99, 064045 (2019), 
arXiv:1812.07865 [gr-qc]. 

Jonathan Blackman, Scott E. Field, Mark A. Scheel, 
Chad R. Galley, Christian D. Ott, Michael Boyle, 
Lawrence E. Kidder, Harald P. Pfeiffer, and Béla Szila- 
gyi, “Numerical relativity waveform surrogate model for 
generically precessing binary black hole mergers,” Phys. 
Rev. D96, 024058 (2017), arXiv:1705.07089 [gr-qc]. 

See Supplemental Material here, for tests of the surro- 


[49] 


[50] 


[51] 


[52] 


[53] 


[54] 


55 


56 


57 


58 


59 


[60] 


[61] 


[62] 


gate models using high spin NR injections. This further 
includes Refs. [57—62]. 

I. M. Romero-Shaw et al., “Bayesian inference for compact 
binary coalescences with bilby: validation and applica- 
tion to the first LIGO-Virgo gravitational-wave transient 
catalogue,” Mon. Not. Roy. Astron. Soc. 499, 3295-3319 
(2020), arXiv:2006.00714 [astro-ph.IM]. 

Vijay Varma, Maximiliano Isi, Sylvia Biscoveanu, Will M. 
Farr, and Salvatore Vitale, “Measuring binary black hole 
orbital-plane spin orientations,” (2021), arXiv:2107.09692 
[astro-ph. HE]. 

Fabio Antonini and Frederic A. Rasio, “Merging black 
hole binaries in galactic nuclei: implications for advanced- 
LIGO detections,” Astrophys. J. 831, 187 (2016), 
arXiv:1606.04889 [astro-ph.HE]. 

G. Monari, B. Famaey, I. Carrillo, T. Pif, M. Steinmetz, 
R. F. G. Wyse, F. Anders, C. Chiappini, and K. Janßen, 
“The escape speed curve of the Galaxy obtained from Gaia 
DR2 implies a heavy Milky Way,” A&A 616, L9 (2018), 
arXiv:1807.04565 [astro-ph.GA]. 

S. Kullback and R. A. Leibler, “On information and suffi- 
ciency,” The Annals of Mathematical Statistics 22, 79-86 
(1951). 

Sizheng Ma, Matthew Giesler, Vijay Varma, Mark A. 
Scheel, and Yanbei Chen, “Universal features of 
gravitational waves emitted by superkick binary black 
hole systems,” Phys. Rev. D 104, 084003 (2021), 
arXiv:2107.04890 [gr-qc]. 

R. Abbott et al. (LIGO Scientific, VIRGO, KAGRA), 
“Tests of General Relativity with GWTC-3,” (2021), 
arXiv:2112.06861 [gr-qc]. 

LIGO Scientific Collaboration and Virgo Collaboration, 
“Gravitational Wave Open Science Center,” https://www. 
gw-openscience.org. 

Michael Boyle et al., “The SXS Collaboration catalog of 
binary black hole simulations,” Class. Quant. Grav. 36, 
195006 (2019), arXiv:1904.04831 [gr-qc]. 

Marissa Walker, Vijay Varma, and Geoffrey Lovelace, 
“Extending numerical relativity surrogate models to near 
extremal spins,” (2021), in preparation. 

Mark A. Scheel, Matthew Giesler, Daniel A. Hemberger, 
Geoffrey Lovelace, Kevin Kuper, Michael Boyle, B. Szila- 
gyi, and Lawrence E. Kidder, “Improved methods for 
simulating nearly extremal binary black holes,” Class. 
Quant. Grav. 32, 105009 (2015), arXiv:1412.1803 [gr-qc]. 
LIGO Scientific Collaboration, Updated Advanced LIGO 
sensitivity design curve, Tech. Rep. (2018) https://dcc. 
ligo.org/LIGO-T1800044/public. 

SXS Collaboration, “The SXS collaboration catalog 
of gravitational waveforms,” http://www.black-holes. 
org/waveforms. 

Katerina Chatziioannou, Geoffrey Lovelace, Michael 
Boyle, Matthew Giesler, Daniel A. Hemberger, Reza 
Katebi, Lawrence E. Kidder, Harald P. Pfeiffer, Mark A. 
Scheel, and Béla Szilagyi, “Measuring the properties 
of nearly extremal black holes with gravitational waves,” 
Phys. Rev. D 98, 044028 (2018), arXiv:1804.03704 [gr-qc]. 


SUPPLEMENTAL MATERIALS 
I. HIGH SPIN NR INJECTIONS 


As GW200129 has a preference for a large x, (see 
Fig. 1), it is important to check the validity of the sur- 
rogate models NRSur7dq4 and NRSur7dq4Remnant in this 
regime. These models were trained on NR simulations 
with x12 < 0.8, but allow extrapolations to x12 = 1 [38]. 
Unfortunately, because NR simulations become expen- 
sive for large spins, very few simulations exist beyond 
X1,2 = 0.8 [57]. Therefore, the surrogate models have 
only been tested against a handful of NR simulations 
outside their training region [30, 38, 58]. While an ex- 
haustive exploration is still prohibitively expensive, we 
conduct two additional tests using high spin NR simula- 
tions from Refs. [54, 59]. We inject these NR waveforms 
(in zero-noise) into a simulated LIGO-Virgo network op- 
erating at design sensitivity [60]. Then, we follow the 
same procedure as described in the main text, and use 
the NRSur7dq4 and NRSur7dq4Remnant models to recover 
the binary parameters and the kick. 

We consider two NR waveforms, SXS:BBH:0178 [59] 
and SXS:BBH:2439 [54], which are publicly available 
through the SXS Catalog [61]. SXS:BBH:0178 corre- 
sponds to an equal mass binary with equal spins that 
are nearly extremal (x1,2 ~ 0.99), and aligned along the 
orbital angular momentum. Because of the symmetries 
(equal masses and spins), the kick for this binary is zero. 
SXS:BBH:2439 also corresponds to a binary with equal 
masses and large spins (x1,2 ~ 0.95), but with most of 
the spin in the orbital plane. In this case, while the spin 
magnitudes are equal, the directions are not. In fact, 


the spins were chosen to lie in the “superkick” config- 
uration [54], with the in-plane spin components being 
anti-parallel to each other; this configuration leads to 
large kicks [4-6]. For both injections, we set dz = 400 
Mpc and jN = Pref = 7/3, where dz is the luminosity 
distance, 6; is the inclination angle between the total 
angular momentum J and line of sight direction N, and 
ret is the reference orbital phase [18]. This leads to large 
signal to noise ratios (SNRs) of 72 and 63, respectively, 
for SXS:BBH:0178 and SXS:BBH:2439, and therefore pro- 
vides a stringent test for the surrogate models. The rest 
of the binary parameters will be shown in figure insets 
below. 


Fig. S1 shows the recovered mass ratio, spins, and kick 
for the two NR injections. In both cases, we find that the 
injected values are well recovered and lie within the 90% 
credible region of the posterior. We note, however, that 
the kick posterior for SXS:BBH:2439, while still capturing 
the true value, peaks below the injected value, which is re- 
lated to the fact that the spin magnitudes also peak below 
the injected value. Similar trends were noted in Ref. [62], 
where it was found that the prior choice of uniform in 
spin magnitude and direction (which we also adopt in this 
work) may not be suitable for binaries with large spins. 
Therefore, an investigation on the impact of prior choices 
on the spin and kick inferences for high spin events may 
be necessary, especially as we approach the high SNRs 
considered in Fig. S1 with future detector improvements. 
In summary, while these tests give us more confidence in 
the constraints placed on GW200129, we recommend a 
more exhaustive study to rule out any potential systemat- 
ics in the surrogate models due to extrapolation to large 
spins, which we leave for future work. 
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Figure $1. Spin extrapolation tests for the NRSur7dq4 and NRSur7dq4Remnant models. Shown are the recovered kick magnitude, 
mass ratio and spin posteriors (measured at tre¢/M = —100) for NR injections. The dark (light) regions represent the 50% 
(90%) credible bounds on joint 2D posteriors, while the diagonal plots show 1D marginalized posteriors. The injected binary 
parameters are shown as black lines, and are also given in the inset text. The left panel corresponds to an equal mass binary 
with nearly extremal (x1,2 ~ 0.99) but aligned spins. The right panel corresponds to an equal mass binary in a superkick 
configuration, with large spins (x1,2 ~ 0.95). In both cases, the injected values for the mass ratio, spins, and the kick are well 
recovered and lie within the 90% credible region of the posterior. As the in-plane spin is zero for the left panel, 61 and ¢2 are 
not meaningful parameters and the offset from the true value is not of concern. 


